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Propagation equations for optical pulses are needed to assist in describing applications in ever 
more extreme situations - including those in metamaterials with linear and nonlinear magnetic 
responses. Here I show how to derive a single first order propagation equation using a minimum 
of approximations and a straightforward "factorization" mathematical scheme. The approach gen- 
erates exact coupled bi-directional equations, after which it is clear that the description can be 
reduced to a single uni-directional first order wave equation by means of a simple "slow evolution" 
approximation, where the optical pulse changes little over the distance of one wavelength. It also 
also allows a direct term-to-term comparison of an exact bi-directional theory with the approximate 
uni-directional theory. 
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I. INTRODUCTION 

In recent years, the propagation of optical pulses un- 
der ever more extreme conditions has been the subject 
of significant attention. This situation has arisen pri- 
marily because of the multitude of applications [l|: e.g. 
where ultrashort pulses are relied upon to act as a kind of 
strobe-lamp to image ultrafast processes 0, 0] , or where 
the electric field profile of a pulse 0, Q is engineered to 
excite specific atomic or molecular responses. Other mo- 
tivations are systems where strong nonlinearity is used 
to construct equally wide-band but also temporally ex- 
tended pulses - i.e. (white light) supercontinua Q-Q ~ 
or even come full circle and use the strong nonlinearity 
to generate sub-structure that is again temporally con- 
fined, as in optical rogue waves @ ; or even the ternpo rally 
and spatially localized filamentation processes [Ifl, 
Further, developments in electromagnetic metamaterials 
p^ - [T^ lead to a requirement for including magnetic dis- 
persion or even magnetic nonlinearity [15| . 

It is clear, therefore, that progress toward shorter pulse 
durations as well as their increasing spectral bandwidths, 
and higher pulse intensities - as well as exotic propaga- 
tion media - are all factors either stretching existing pulse 
propagation models to their limits, or breaking them. In 
such regimes, we need to be sure that our numerical mod- 
els still work, and have a clear idea of what has been 
neglected, and what the side-effects of those approxima- 
tions are. Most existing pulse propagation models make 
sequential approximations that can have unforeseen side 
effects. In contrast, in this article, I show how a straight- 
forward and relatively simple derivation allows a side- 
by-side comparison of exact and approximate propaga- 
tion equations, whilst still providing the numerical and 
analytical convenience of a first-order wave equation. 

The analysis of optical pulse propagation traditionally 
involves describing a pulse in terms of a complex field en- 
velope, while neglecting the underlying rapid oscillations 
at its carrier frequency. The resulting "slowly varying 
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envelope approximation" (SVEA) (see e.g. which 
reduces second order differential equations to first order, 
is valid when the envelope encompasses many cycles of 
the optical field and varies slowly. Starting with the sec- 
ond order wave equation, other auxiliary assumptions are 
required to get the final result of a first-order wave equa- 
tion: the introduction of a co-moving frame, and the 
neglect of usually negligible second order spatial deriva- 
tives. Although it is now easily possible to choose to solve 
Maxwell's equations numerically instead (see e.g. (itI - 
[2l|). the approach lacks the intuitive picture of a pulse 
"envelope" , and tends to be computationally demanding. 

Many attempts have been made to generalize the 
SVEA style of derivation, and perhaps the most no- 
table of these was that of Brabec and Krausz [iQl. By 
slightly relaxing one assumption, they derived corrections 
to the SVEA, which they included in their "slowly evolv- 
ing wave approximation" (SEWA). This enabled the few- 
cycle regime to be modeled with improved accuracy, and 
the SEWA has subsequently been applied in different sit- 
uations, including ultrashort IR laser pulses in fused sil- 
ica 123^ [2J|, the filamentation of ultra-short laser pulses in 
air '24] , and even in micro-structured optical fibres [2^ . 
Later, Porras [2^ proposed a slightly different "slowly 
evolving envelope approximation" (SEEA) that included 
corrections for the transverse behavior of the field; and 
Kinsler and New [l^ took the process as far as it would 
go with their "generalized few-cycle envelope approxima- 
tion" (GFEA). Although the wave equation generated by 
the GFEA was generally too complicated for practical 
use, its derivation exposes one important point: extend- 
ing SVEA style derivations into wide-band situations ex- 
poses the user to a number of poorly controlled side ef- 
fects ^2§\ . Many other styles of dervivation also exist (see 
e.g. [2t^31.] ). but most use similar approximations, and 
apply them sequentially. 

Here I will show that an alternative "factorization" 
style of derivation we can achieve the simplicity of a first- 
order wave equation for optical pulse propagation, but 
avoid the unpleasant side-effects of the traditional ap- 
proach. Early but rather limited examples are by Shen 
|16l | , Blow and Wood ^ , and perhaps Husakou and Her- 
rmann [ssj ; more recently (and more rigorously) we have 
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Ferrando et al. [3J| and Genty et al. (35[. The math- 
ematical basis of the factorization shown in this article 
relies on Ferrando et al. [s^l, but here I make a point 
of generating wave equations incorporating most optical 
effects - both electric and magnetic dispersion, diffrac- 
tion, second and third order nonlinearity, angle depen- 
dent refractive indices, and so on. In particular, prior 
to any approximations being applied, there is an (explic- 
itly bi-directional) stage where two counter-propagating 
wave equations are coupled together. This provides us 
with an important insight: that a simple "slow evolu- 
tion" approximation is all that is needed to obtain a uni- 
directional first order wave equation, irrespective of the 
origin of the coupling. 

In this article I give a description of a modern approach 
to optical pulse propagation applicable to most situations 
that occur in nonlinear optics. This is a regime where 
we want to model the most general situations possible, 
while avoiding having to do a full numerical simulation 
of Maxwell's equations. The treatment here is intended 
to be straightforward enough for the student, whilst also 
being comprehensive enough so that both novice and spe- 
cialist can really understand the nature and limitations of 
this and other pulse propagation models. Starting with 
a general form of the second order wave equation in sec- 
tion |TT1 I follow with discussion the important role of the 
choice of propagation direction in section IIII[ which in 
nonlinear optics is usually in space and not in time. In 
section llVi I introduce the method of factorization that 
allows us to construct an explicitly bi-directional model, 
and which is then reduced to the uni-directional limit in 
section |Vl where nonlinear pulse propagation is typically 
applied. Section IVTl discusses typical modifications that 
can be applied to the equations given in sections IIVI and 
IVI in order to and simplify them appropriately and com- 
pare them to existing models; whilst section I VIII gives 
specifc examples for the common cases of propagation 
media with second and third order nonlinearities. The 
article is then summarized in section rVIIII 



II. SECOND ORDER WAVE EQUATION 

Most optical pulse problems consider a uniform and 
source free dielectric medium. In such cases a good 
starting point is the second order wave equation for 
the electric field, which results from the substitution 
of the \7 X H = dtD + J Maxwell's equation into the 
\/ X E = —dfB one (see e.g. [Ill), although here I also 
allow for free currents J. Magnetic effects can also be 
incorporated - easily so in the case of linear magnetic 
dispersion, but also it is possible to retain a term for 
more general magnetic effects. However, cases where ei- 
ther the permittivity e{uj) or the permeability /i(w) are 
negative are not excluded. 

A sufficiently general model of the dielectric response 



in the time domain is 

D{r,t)=e{t)^E{r,t) (1) 
= eofi [r, t)^E{r,t) + eoPe{E,f,t), (2) 



where the scalar e^, contains the linear response of the 
material that is both isotropic^ and lossless (or gain- less); 
since here it is a time-response function, it is convolved 
with the electric field E. Note that the field vectors E, D, 
and indeed the material parameter are all functions 
of time t and space r = (x.y, z); the polarization is 
a function of time i, space r, and field E. The following 
derivation also allows for magneto-electric polarizations, 
i.e. those where also depends on H, although I do 
not explicitly include such a dependence in the notation. 
Similarly, the magnetization response is 



Bif,t)^tiit)*Hif,t) (3) 
= Hol^Lir, t) * H{f, t) + fioMf^iH, r, t), (4) 



where the scalar fi^ contains the linear response of the 
material that is both isotropic and lossless (or gain- less). 
Note that H,B and fiL are all functions of time t and 
space r = {x,y,z); the magnetization is a function 
of time t, space r, and field H. The following deriva- 
tion also allows for magneto-electric magnetizations, i.e. 
those where also depends on E, although I do not 
explicitly include such a dependence in the notation. 

Since here I have chosen to incorporate the "simple" 
linear responses of the propagation medium in and 
/ii, the remaining parts P^, will usually be in part 
electric and magnetic field dependent, and incorporate 
effects such as birefringence, angle dependence, and non- 
linearity; it should also incorporate any loss 37[. For 
example, P^ might contain a scalar nonlinearity such as 
third order Kerr nonlinearity with P„/ oc {E ■ E)E, or 
a (vector) second order nonlinearity. Note that it is not 
always necessary or desirable to include all the simple lin- 
ear responses in and /ii, some may be left in P^, M^; 
as will be discussed later. Alternatively, and in accor- 
dance with [s^l we could choose to pick cl and such 
that el/zl is real, rather than each being real valued on 
its own. However, this would alter the handling of the J, 
Pf, and terms. 



Defining V = {dx,dy, dz) and da = d/da, eoMo = 
and current density J, we can write the exact second 



The isotropy of (and later of /j^) is both important and uscfuL 
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order wave equation as 



c^V X V X E{t) = -d'f 



^lL{t)^eL{t)*E{t) 



y.otJ'L{t)*dtJ - dt 



V X Af^(t) 



(5) 



Here I have suppressed the space coordinates and electric 
field dependence for notational simplicity. The (usual) 
next step is to replace V xV x E above with the identity 
VV-i^-V^-E, where as usual S/"^ = d^+dl + d^. Initially 

this might look over-complicated, since VV • E adds in 
some extra terms (e.g. a d'^E^) which are then canceled 
by the same term from V^i?. However, since the field 
divergence is an important Maxwell's equation, splitting 
the double curl operation in this way turns out to be 
advantageous. 

For the case of a free charge density p, and with the 
same separation of the material response as used above, 
Maxwell's equations tell us that 



V-L' = /9 = eoV- eLirE + P, 



(6) 



= eoeL *V-E + eo [Vel] • -kE + eoV ■ P, 

(7) 



(8) 



so for an isotropic e^, we can use Ve^ = 0; note that 
isotropy also implies field-independence. The frequency 
domain changes convolutions into products, so that we 
have 

eoeL(^)V-^(cj) = p(w)-V-Pe(w) (9) 
p(w) V-Pe(c^) 



V • E{uj) 



(10) 



Note that the left-hand side (LHS) of this equation (i.e. 
V • E) seems to be potentially large, since it consists of 
field derivatives. However, the divergence condition re- 
veals that with no free charge it is simply V-Pc/ei,, which 
merely is of the order of the nonlinearity or anisotropy 
of e; both of which are small in typical systems. Since 
V[V • E] is typically much smaller than V^E, it can rea- 
sonably be considered as a correction to a propagation 
dominated by V^E. 

As a result, we find that the replacement of V x V x 
E by —V^E + VV • E not only achieves this valuable 
minimization, but it also reduces the remaining spatial 
derivatives to the simple \/^E. The side effect is that 
we now need to compute VV • P^, which may well be a 
complicated function of E; it also gives rise to phenomena 
such as nonlinear diffraction term (see e.g. (38|). 



The second order wave equation is best written in the 
frequency domain, because of the need to divide the di- 
vergence term by the frequency dependent el] and so is 

-c^V^Eiu) ^ Lu^d^eL{uj)nLiuj)E{uj) + Lo^fiL{i^)P,{uj) 

+ iw/XL(a;) /(w) + ^— V x 
Co 

„2 



eL(w) 



V-P,iuj)- 



(11) 



For plane polarized pulses, a scalar version allowing for 
just one of the linear polarization components is suffi- 
cient. However for materials that couple the horizon- 
tal and perpendicular polarizations together, such as the 
X*-^' interaction relied on by optical parametric amplifiers 
(OPA) or oscillators (see e.g. [l^), we could write one 
equation for each polarization, and then find that they 
were coupled together by the nonlinearity. 

The wave equation in eqn. (jlip contains both current 
J and charge density p terms, which are usually interde- 
pendent. These terms are not often important in pulse 
propagation, so I do not discuss their modeling; appro- 
priate treatments can be seen in the literature on optical 
filamentation (see e.g. (40l|V 



III. PROPAGATION DIRECTION 

In this article I will not be considering strong reflec- 
tions from material modulations or interfaces. Never- 
theless, considering simple reflections is an excellent way 
of clarifying some important issues that arise when we 
choose whether to propagate pulses forward in time, or 
forward in space. 

Temporal propagation is the usual choice in flnite 
difference time domain (FDTD) modeling of Maxwell's 
equations 0, HI]) where fields E{x,y,z),H{x,y,z) are 
stepped forward in time t; exitations of the fleld (i.e. op- 
tical pulses) then evolves backward or forwards in the 
space coordinates {x,y,z). We therefore set up initial 
conditions covering each point in space at a chosen ini- 
tial time ti] likewise we read out our final state for each 
point in space at a chosen final time t/, as shown on fig. 
[T] This choice requires a time-response treatment of dis- 
persion, perhaps involving convolutions, however, as also 
shown by fig. [U it provides natural reflections. 

Spatial propagation is the usual choice in nonlin- 
ear optics and optical pulse propagation, where flelds 
E{t, X, y), H{t, X, y) are stepped forward in a chosen spa- 
tial direction (z); exitations of the fleld (i.e. optical 
pulses) then evolves backward or forwards in time and 
space coordinates {t,x,y). We therefore set up initial 
conditions covering each point in time at a chosen point 
in space Zi] likewise we read out our final state for each 
point in time at a chosen point in space zy, as shown 
on fig. [2j Comparison of figs. [1] and [2] also show that 
to be correctly modeled, an ordinary refiection from the 
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FIG. 1: An ordinary reflection at an interface between media 
witli permittivities ei and £2, in a t-propagated picture. An 
incoming pulse propagates forward (in t) and evolves forward 
(in z) until it reaches an interface, whereupon it splits into a 
transmitted pulse and a normal reflected pulse; the reflected 
pulse then evolves backward in space as both transmitted and 
reflected pulses continue to propagate forward in time. 



interface back to our initial point must be included in 
our initial conditions. Unfortunately, we will usually not 
know the properties of this reflection in advance, so we 
will not include it in the initial conditions. As a result, 
our solution of Maxwell's equations at the interface cre- 
ates the mirror image pulse that is needed to exactly 
cancel out the ordinary reflection. Next, since we have 
chosen to propagate solely toward larger z, this mirror 
image '"reverse reflection" pulse now evolves forward in 
space z but backwards in time as shown on fig. [51 

We see, therefore, that if we want to take advantage 
of the benefits of spatial propagation, notably the efi- 
cient handling of dispersion, we will also not want to be 
modeling systems containing significant reflections. In- 
deed, this issue motivated the time-propagated model of 
Scalora et al. [13, S, HI] , which are based on the second 
order wave equation; however that approach suffers some 
of the same drawbacks as other tradition pulse propa- 
gation techniques. To handle a temporally propagated 
model based on a second order wave equation, it is best 
to use that for the displacement field D rather than for E; 
since time derivatives of D appear directly in Maxwell's 
equations, whereas those for E are complicated by the 
material response. 



A. Spatial propagation 

The first step to achieving a first order wave equation 
containing the necessary physics but without unnecessar- 
ily complex approximations is to reorganize the wave eqn. 
([5]) to emphasize contributions that by themselves can 
freely propagate forward and backward without interact- 
ing. To do this I choose a specific propagation direction 



FIG. 2: A reflection at an interface between media with per- 
mittivities ei and e2 , in a ^-propagated picture. An incoming 
pulse propagates forward (in 2) and evolves forward (in t) 
until it reaches an interface, whereupon it splits into a trans- 
mitted pulse and its reverse reflection; the reverse reflected 
pulse then evolves backward in time as both transmitted and 
reflected pulses continue to propagate forward in space. A re- 
verse reflection is the means by which a spatially propagated 
system represents a pulse propagating backwards in z, so that 
it cancels the ordinary reflection missing from the initial con- 
ditions. 



(e.g. along the z-axis), and then denote the orthogonal 
components (i.e. along x and y) as transverse behaviour; 
many situations are also cylindrically symmetric, allow- 
ing simplification of the two transverse dimensions x,y 
into a single radial coordinate r. I therefore rearrange 
eqn. pT|) into 



— ikQ^i^c^^ J{uj) — ikocV X 

' -[v.p.(.)-^(") 



V 



(12) 



where fcg = lo'^/c^ and P'^iuj) 
uj'^eoijLoeL{oj)iiL{'^)', k^ — uP' jc?. Here all the simple 
linear response (e.g. the isotropic refractive index and 
dispersion) has been moved to the LHS as a (possibly) 
frequency dependent propagation wave vector; the resid- 
ual responses (i.e. and M^) contain any non-oj de- 
pendence, angle dependent terms, nonlinearity or spatial 
variation. Note that defining /3(aj) is a matter of choice^ 
in some cases we may find it convenient to define it to 
be frequency independent; in others we might (e.g.) even 
decide to retain some angle dependence, perhaps even to 
the point of generating a spherical "in-out" bi-directional 
model, rather than a linear forward-backward one. 
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IV. FACTORIZATION 

I now factorize the wave equation, a p rocess which, 
while used in optics for some time 1321 has only re- 
cently been used to its full potential [3J, [H, . Fac- 
torization neatly avoids almost all of the approximations 
necessary in the standard approach and its extensions 
[19I [27L [31I Isgj (etc) - which are in fact much more com- 
plicated than they first appear, as has been shown by 
detailed analysis [23, H^. A major advantage of fac- 
torization is that we can directly compare the exact bi- 
directional and approximate uni-directional theories term 
for term, whereas in other approaches the backward parts 
simply vanish and are not directly available for com- 
parison. Perhaps the clearest recent description of the 
approximations made in a standard (non-factorization) 
derivation of a uni-directional wave equation is by Berge 
and Skupin [i^. That work discussed the filamentation 
resulting from nonlinear self-focusing effects, so that they 
incorporated the role of the longitudinal field components 
and included a model for a plasma connecting the J and 
p contributions - here I retain these terms, but the reader 
is referred to Berge and Skupin [i^ for a specific model. 

Factorization takes its name from the fact that the 
LHS of eqn. (jl2p is a simple sum of squares which might 
be factorized, indeed this is what was done in 1989 in a 
somewhat ad hoc fashion by Blow and Wood [s^ . Since 
the factors are just (9^ =Fz/?), each by itself looks like a for- 
ward (or backward) directed wave equation. A rigorous 
factorization procedure [sS |46| , of which some basics are 
given in appendix |X1 allows us to define a pair of counter- 
propagating Greens functions, and so divide the second 
order wave equation into a bi-directional pair of coupled 
first order wave equations. That these factorized equa- 
tions are equivalent to the original second order wave 
equation is proven by taking their sum and differences, 
then substituting one into another with the assistance 
of a derivative with respect to explained in Ref. 

[s^l^. Further, even in the approximate uni-directional 
limit, the factorised wave equations have been shown by 
Genty et al. [35| to give a stunning level of agreement 
with pseudospectral spatial domain (PSSD) [2l| Maxwell 
equations simulations. 

Before proceeding, it is worth reiterating an important 
point - the choice of £1.(0;) and /iL(w), and therefore of 
/3(w) in eqn. ([T^ . defines the specific Greens functions 
used; it therefore also defines the underlying basis on 
which we will then propagate the electric field E. 

As an aside, the interested reader may wish to examine 
the mathematical "wave-splitting" work of Weston and 
others (see e.g. [13]), although it does not consider resid- 
ual terms, and (at least initially) was primarily concerned 
only with reflections and scattering. This was based on 
that from the earlier work of Beezley and Krueger [i^ 



who applied wave-splitting concepts to optics. Other sim- 
ilar work is the one-way wave equation of Leviandier Boj ] , 
and other directional schemes have been suggested by 
Kinsler et al. [13] and Kolesik et al. [Fl'l . It is also inter- 
esting to compare and contrast the factorization scheme 
used here with beam propagation methods (BPM, e.g., 
Refs. I^qI, ^S^j). For example, the treatment of Van Roey 
et al. [29] also begins using Greens' functions which de- 
flne a chosen reference propagation. Thus, whilst those 
BPM methods might (in principle) be developed in a way 
which matches the benefits of the factorization scheme I 
present here, to my knowledge no such implementation 
has been published. 



A. Bi-directional wave equations 

A pair of bi-directional wave equations suggests simi- 
larly bi-directional fields, so I split the electric field into 
forward {E~^) and backward {E~) directed parts, with 
E = E^ + E^ . In the following equations I have rein- 
stated the E argument of (and H of M^) to emphasise 
that they depend on the total field; an important point 
since we see that P^, M^, diffraction, and other terms 
drive both forward and backward equations equally. 

Using the procedure summarized in Appendix |21 the 
second order wave equation in eqn. (1121) can be converted 
into a pair of coupled bi-directional first order wave equa- 
tions for the directed fields E^ . They arc 



dzE^{ijj) = ±iP{uj)E^{uj) ± 



E+ {lu) + E- {uj) 



± 



ikl{uj)iiL 
2/3(a;)c 



± 



2/3(cj)ei(w) 



P,{E+ +E-,E,,u) 
2/3(w) 

V-F,(w)- 



p(w) 



(13) 



See section IV. B of this reference. 



Since ko = uj/c, such factors convert to a (scaled) time 
derivative when these frequency domain equations are 
transformed into the time domain. 



B. Propagation, evolution, and directed fields 

Note that since our solutions of the wave equations 
enforce propagation toward larger z, the fields E^{t) are 
directed forwards and backward in time; these fields then 
evolve forwards and/or backward in time as z increases. 
Note that I use this terminology (propagated, directed, 
evolved) throughout this article to mean these three spe- 
cific and distinct concepts. 

When examining the wave equation eqn. (|13l) which 
evolves the directed fields E^ as they propagate forward 
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in z, we see that the right-hand side (RHS) has two types 
of terms: which I label the "underlying" and "residual" 
parts [37| . 

Underlying evolution is that given by ±il3{uj)E^ term, 
and is determined by our chosen eL('^) and ijll(uj). By 
itself, it would describe a plane- wave like evolution where 
the field oscillations would move forward or backward 
(— ) in time across E^{t). This is analogous to the choice 
of reference when constructing directional fiel ds [ 50l | , or 
the refractive index term Uq used in the BPM [29| . 

Residual evolution accounts for the discrepancy be- 
tween the true evolution and the underlying evolution, 
and is every part of the material response not included 
in eL{^) or i.e. it is all the remaining terms on 

the RHS of eqn. (fT3|) . These typically include any non- 
linear polarization, angle dependent linear terms, and the 
transverse effects; they are analogous to the correction 
terms used in directional fields models, or the refractive 
index perturbation An^ used in BPM [29]. In the lan- 
guage used by Ferrando et al. Q, these residuals are 
"source" terms. Although we might hope they will be 
a weak perturbation, so that we could make the (desir- 
able) uni-directional approximation discussed later, the 
factorization procedure is valid for any strength. 



C. Underlying evolution: choice of /? and the 
resulting 

I now examine how the choice of /3 affects the relative 
sizes of the forward and backward directed and E' . 
To do this I consider the simple example of a medium for 
which the field is known to propagate with wave vector fc; 
but for demonstration purposes we choose an underlying 
evolution determined by a wave vector /3 that is differ- 
ent from k. For example, for a linear isotropic medium 
we could exactly define = -I- A^; but in general 
we would just have some residual (source) term Q. This 
means that our definitions of forward and backward di- 
rected fields do not exactly correspond to what the wave 
equation will actually evolve forward and backward as 
we propagate toward larger z. 

The second order wave equation is (9^ -f /3^)i? = — Q, 
which in the linear case has Q = S'^E, so that {d^ + 
k'^)E = 0. The factorization in terms of /3 is then 



^ 213 



(14) 



Now if we select the case where our field E only 
evolves forward, we know that E = i5o exp[«A:2;]. Conse- 
quently E'^ must have matching oscillations: i.e. E^ — 
E^ exp[ikz], even though E~ is directed backward. Sub- 
stituting these into eqn. gives 



E. 



(3 



P + k 



(15) 



which specifies how much E we need to combine with 
E^ so that our pulse evolves forward; since the E~ will 



be dragged forward by its coupling to E'^ . This interde- 
pendent E^ behaviour is generic - no matter what the 
origin of the discrepancy between /? and the true evolu- 
tion of the field (i.e. the residual or source terms such 
as mismatched dispersion, nonlinearity, diffraction, etc): 
some non-zero backward directed field E^ must exist but 
still evolve forwards with E^ . Analogous behaviour can 
be seen in the directional fields approach of Kinsler et al. 

Usually we hope that this residual E contribution is 
small enough so that it can be neglected. If we assume 
E- ~ 0, then we find that k ~ (3 + A2/2^, which is 
just the expansion of A: = (/3^ -I- A^)-'^/^ to first order in 
A^//3^. Following this, we find that eqn. dTSl) then says 
that E^ ~ (A^/4/3^)_E(|, which has come fuU circle and 
provided us with the scale on which E~ can be considered 
negligible. Outside the restricted (linear) case where we 
know A^ , the true wave vector k might be difficult to de- 
termine, and in nonlinear propagation may even change 
as the pulse propagates. 

There is a further important point to notice: if we 
choose j3 — l3{uj) with a frequency dependence, then we 
see that the source-like terms (e.g. diffraction, polariza- 
tion, etc; or A^ in eqn. (|T4l) ') inherit that dispersion. This 
means that even if we started with polarization model 
with instantaneous nonlinearity, our factorized equations 
no longer have instantaneous nonlinear terms, as they 
have become "anti-dispersed" by the factor of /3(a;)~^; 
as indeed have the other residual terms. This matches 
exactly what happens in the directional fields approach 
of Kinsler et al. [50|, where choosing a dispersive refer- 
ence has an equivalent effect on the correction terms. 



V. UNI-DIRECTIONAL WAVE EQUATIONS 

Making only a single well defined type of approxima- 
tion I can now reduce the exact coupled bi-directional 
evolution of E± down to a single uni-directional first or- 
der wave equation. I do not require a moving frame, a 
smooth envelope, or to assume inconvenient second or- 
der derivatives are somehow negligible: all these are fre- 
quently required in standard treatments, and even exten- 
sions use them lli l26i, 113, [S il, 39]. The approxima- 
tion is that the residual terms are weak compared to the 
(underlying) ±i/3E term - e.g. weak nonlinearity, angle 
dependence, and diffraction. This enables me to assert 
that if I start with E^ = 0, then E~ will remain negligi- 
ble - see my estimate in subsection II V CI In this context, 
"weak" means that no significant change in the backward 
field is generated in a distance shorter than one wave pe- 
riod ("slow evolution"); and that small effects do not 
build up gradually over propagation distances of many 
wavelengths ("no accumulation"). 

Slow evolution is where the size of the residual terms 
is much smaller than that of the underlying linear evolu- 
tion - i.e. smaller than /SE. This allows us to write down 
straightforward inequalities which need to be satisfied. It 
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is important to note the close relationship between these 
and a good choice of /?, as discussed in subsection II V CI If 
^ is not a good enough match, there always be significant 
contributions from both forward and backward directed 
fields; and even if nothing ends up evolving backwards, 
an ignored backward directed field will result in miscalcu- 
lated nonlinear effects, since the total field E = 
will be different to the assumed value of . 

No accumulation occurs when the evolution of any 
backward directed field E^ is dominated by its coupling 
via the residual terms to the forward directed field E^; 
and not by its preferred underlying backward evolution. 
No accumulation means that forward evolving field com- 
ponents do not couple to field components that evolve 
backward; this the typical behaviour since the phase mis- 
match between forward evolving and backward evolving 
components is ~ 2/3; in essence it is comparable to the 
common rotating wave approximation (RWA). This rapid 
relative oscillation means that backward evolving compo- 
nents never accumulate, as each new addition will be out 
of phase with the previous one; it is not quite a "no reflec- 
tion" approximation, but one that asserts that the many 
micro-refiections will not combine to produce something 
significant. An estimate of the conditions required to 
break this approximation are given in Appendix |B] gen- 
erally speaking this is a much more robust approximation 
than the slow evolution one. Of course, periodic spatial 
modulation of the medium gives periodic residual terms, 
and these can be engineered to force phase matching. In 
most contexts this would be a periodicity based on a rela- 
tively small phase mismatch (see e.g. quasi phase match- 
ing in Boyd [1^); but might even go as far as matching 
the backward wave (see e.g. [ssjV 

It is also important to note that the same small size of 
perturbation from the residual terms can accumulate on 
the forward evolving field components (or, indeed, the 
backward perturbation on the backward evolving field 
components). Although the magnitude of the residual 
terms acting on the forward and backward field evolution 
are identical, forward evolving components of the residu- 
als can accumulate on the forward evolving field because 
they are phase matched; whereas backward residuals are 
not, and rapidly average to zero. 



A. Polarization and Magnetization 

To see most clearly how different optical effects satisfy 
this slow evolution criteria, I will split the total polariza- 
tion Pe into pieces: 

fiL{t)*PeiE,f,t)^ME,t)* E{r, t) + Ve{E, r, t) 
= H{E.t)^E{r,t) 
+ (j)N{E,t)^E{r,t) 
+ VLiE,r,t) + VNiE,r,t). (16) 



The part which is scalar in nature is represented by 0e, 
it might contain linear parts and time response {(t>L)', 
but can also be a function of transverse wave vector 
(i.e. be angle dependent), or contain nonlinear contribu- 
tions 4>N such as the third order Kerr nonlinearity with 
(j)NE oc {E ■ E)E. The vector part would typically be 
e.g. a second order nonlinearity, which couples the or- 
dinary and extra-ordinary field polarizations. Note that 
this description of the material parameters does not re- 
strict allowed values of e in any way; they can include 
any order of nonlinearity. 

The same can be done for M^, the non- isotropic and 
nonlinear (i.e. the non-/ii) part of the magnetization. 
However, the calculations will all follow the same basic 
pattern that they do for , albeit somewhat complicated 
by the curl operation. Since magnetic nonlinearity is 
rarely present when considering optical propagation, I 
leave detailed assessment of such effects to later work. 



B. Residual terms and slow evolution 

Now I will treat each possible residual term in order, 
where the oppositely directed field is negligible: i.e., for 
E^, we have that E^ ~ 0. where the scalar con- 
tains the linear response of the material that is both 
isotropic and lossless (or gain- less); since here it is a time- 
response function, it is convolved with the electric field 
E. Note that the field vectors E, D, and indeed the ma- 
terial parameter are all functions of time t and space 
f= {x,y,z)] the polarization P^ and its components ^e, 

are a functions of time t, space f, and the field E. 

Below I will refer to field components Ei, where E = 
{Ex, Ey, Ez) and i G {x,y,z}; also to wave vector com- 
ponents ki from k — {kx,ky,kz), with fc^ — k^ + ky. 
However, note that in the constraints below, that fcx is 
also used as a substitute symbol to represent any one of 
kx, ky, or 

Firstly, we have the diffraction term V\E, which is lin- 
ear. For j,j G {x,?/}, and in transverse wave vector 
space, the criteria is 



^k- 


Et+E- 


1/2/3 _^ tk^\ 




/2/3 




il3\Et 





This is just the criterion already given in [5J|, and is 
identical to the standard paraxial criteria. This diffrac- 
tion constraint applies only to the transverse behaviour 
of the pulse, it does not constrain the pulse's intensity, 
temporal bandwidth, or field profile in any way. 

Secondly, scalar polarization terms 4>e, which can be ei- 
ther linear (0l) or nonlinear ((f>]\[{E)). These might en- 
code e.g. some of the dispersion, birefringence, or per- 
haps an angle-dependent refractive index; if nonlinear 
they might arise from e.g. a third-order nonlinearity. 
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Such terms give us the criterion 









1/2/3 _ 




/2/3 


iP\Ef\ 


\Ef 





2/32 



< 1. 



(18) 



In the linear case, (f> = (p^ is independent of E, so only the 
material parameters are constrained, the pulse properties 
play no role. In the nonlinear case, e.g. for a third- 
order nonlinearity, as already treated in (sBI . l45j , we have 
(f> = ~ x^^M^'^l'^- Thus the nonlinear criteria makes 
demands on the peak intensity of the pulse - but does not 
apply smoothness assumptions or bandwidth restrictions. 

Thirdly, linear and nonlinear terms from V^. These will 
have a criterion broadly the same as the scalar cases 
in eqn. (IT5|) . but with V^^ replacing (p^E. Thus for 
i (z {x, y, z}, we can write down constraints for each com- 
ponent of the vector Vg, which are 



IK. 



Ef 



(19) 



In the linear case, K = Vl, and since Vl and E have 
some linear relationship, this criterion only constrains the 
material parameters contained in Vl, not the pulse. In 
the nonlinear case, K = Vn, the same holds except just 
as for scalar nonlinear terms, the peak pulse intensity is 
restricted; e.g. for a x^^^ medium, \Vn\ x^^^l^l- 

However, one complication of the vector cases is that 
a field consisting of only one field polarization compo- 
nent (e.g. E^) may induce a driving in the orthogonal 
(and initially zero) components (e.g. E^). Hence both 
Ey fields will be driven with the same strength, so that 
it is far from obvious that we can set E~ to zero, but 
still keep the E^ without being inconsistency. However, 
as described above, it is the phase matching which en- 
sures that forward residuals accumulate, whilst the non- 
matched backward residuals are subject to the RWA, and 
become negligible: hence we can still rely on eqn. (|19p. 
albeit under caution. 

Fourthly, we have the divergence term VV-Pe. Often this 
term is considered negligible, and discarded even before 
writing down the second order wave equation; neverthe- 
less we should test it. Here we consider just scalar linear 
or nonlinear terms 4>, but the arguments can be adapted 
to vector terms as done above; in any case the results are 
comparable. For i,j £ {x,y, z}, we have 



kikj 
2^^ 



+ Pr|/2/3<z;3|S=' 
\E+ + E-\ 



(20) 



There are four distinct cases to consider here, but only 
two resulting criteria. First, if i G {x,y}, then whether 



j e {x, y} or j — z we find that 



2/32 



< 1 



(21) 



since [E'^j/jiJil ~ k±/l3; this we see that this is a combi- 
nation of both the diffraction and nonlinear criteria, and 
is thus easily satisfied. For the second, where i = z, all 
the wave vector contributions cancel, leaving simply 



< 1. 



(22) 



It is thus directly comparable to the scalar nonlinear cri- 
teria above, and equally likely to be satisfied; the com- 
parable vector criteria are fc^Vi/2/32 ^ Ei and Vi <^ Ei. 

Fifthly, we must consider the charge density p and charge 
current J. These criteria are simple to write down, but 
whether they are satisfied will depend on the initial con- 
ditions and the response of how these are modeled to the 
propagating pulse. This is something that may need to be 
checked during simulation or solution of the pulse propa- 
gation, and not assumed beforehand, although Berge and 
Skupin [4^ discuss the issues in the context of optical 
beam filamentation. The charge and current constraints 
are 



2/32eo H 
koHL \Ji 



2/32c 



< \E, 



« \E, 



(23) 
(24) 



Sixthly, a constraint on the non-^u^ magnetization AI^ 
can also be written down, although (as already discussed) 
I leave the details for later work. It is 



kpc 



V X M„ 



« \E^ 



(25) 



Here the curl operator might often be expected to return 



a value of order /3, so with fco 
\EA. 



13 we have c\M^\/2 < 



To summarize, the diffraction criterion asserts the beam 
must be sufficiently paraxial, the linear criteria asserts 
the material must have weak dispersion, and the non- 
linear criteria assert the nonlinear effect must be weak. 
Paraxiality is determined by our experimental conditions, 
and can thus be guaranteed if desired, and for most opti- 
cal materials, the dispersion is sufficiently weak - except 
perhaps in the vicinity of resonances or band gaps. Weak 
nonlinearity is invariably guaranteed by material damage 
thresholds, since the material suffers damage long before 
nonlinear effects become strong - nevertheless, the effects 
of such strong nonlinearities on uni-directional approxi- 
mations have been analytically and numerically studied 
[4^. Finally, it is worth noting that each criterion is in- 
dependent of the others, so each effect can be tested for 
separately. 
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C. Uni-directional equation for E'^ 

In the case where all of the wavelength-scale slow- 
evolution criteria listed above hold, we can be sure that 
the backward directed field E~ is negligible, and if the 
no-accumulation condition also holds, then neither will 
there be any backward evolving contributions to the field. 
Consequently, we can be sure that an initially negligi- 
ble E~ remains so, and again with ko — lo/c, the bi- 
directional eqn. ((T5)) simplifies to 



2/3(w) 



2/3(a;)c 



V X M^{H^ 



2p{u)tL{uj) 



V 



V-f5(^+, 



£0 



(26) 



Here now the polarization , diffraction, and divergence 
are solely dependent on the forward directed field [E^). 
Likewise the magnetization term should be consid- 
ered as being solely dependent on the forward directed 
field {H'^) - although we will need to estimate the value 
of using the known electric field i?"*" . Since we are in 
a slow evolution approximation, a good estimate for the 
components of will simply be those of i?+ scaled by 
eo(eL/ML)^^^c; so that (e.g.) depends on Also, 

the V X will be dominated by the z dependence of 
its X and y components, so that it will typically generate 
factors of order f3\M^\. 

Although I have included magnetic effects in the 
derivation of eqn. (I^Bl) . I do not consider specific cases 
in detail, as has been done for plane-polarized light in 
e.g. [i35-57]. The derivations in those articles are "tradi- 
tional" in the sense that each consists of multiple interim 
stages at which an additional approximation is applied; 
it is instructive to compare those derivations with mine. 
In particular, e.g. all apply bandwidth limitations, and 
discard various high-order derivative terms that are not 
specific to their choice of propagation medium. Although 
Scalora et al. [55| is the least aggressive in this respect, 
it does not allow for magnetic nonlinearity. 



VI. MODIFICATIONS 

Let us now consider some of the strategies used in 
other approaches, some of which were required in order 
to get approximations that eventually achieved a suffi- 
ciently simple evolution equation. In p articular, the var- 
ious envelope equations (e.g. [l^ [UHO, HH, and even 
[sgI [131) all use co- moving frames and/or envelopes as a 
preparation for discarding inconvenient derivatives: here 



such steps are optional extras. In this factorization ap- 
proach shown here, none of these were required, but they 
nevertheless may be useful. Examples are as follows: 

1. A co-moving frame can now be added, using t' = t— 
z/vf. This is a simple linear process that causes no 
extra complications; the leading RHS if3E^ term 
is replaced by i(/3 =F kf )E^, for frame speed Vf = 
LUi/kf. Note that setting (3 — kf will freeze the 
phase velocity of a pulse centred at lui, not the 
group velocity. 

2. The field can be split up into pieces localized at cer- 
tain frequencies, as done in descriptions of OPAs 
or Raman combs (as in e.g. [13, [H, [111). The 
wave equation can then be separated into one 
equation for each piece, coupled by the appropri- 
ate frequency-matched polarization terms (see e.g. 

3. A carrier- envelope description of the field is 
not required, but can easily be implemented 
with the usual prescription of [s^, [61 1 E{t) = 
A{t) exp[i{ujit — kiz)] -f A* {t) exp[—i{ujit — kiz)] 
defining the envelope A{t) with respect to carrier 
frequency uji and wave vector ki ; this also provides 
a built in a co-moving frame Vf — wi/fci. Multi- 
ple envelopes centred at different carrier frequencies 
and wave vectors (wi, ki) can also be used |39l. [60|. 

4. Bandwidth restrictions might be added (see below), 
either to ensure a smooth envelope or to simplify 
the wave equations; in addition they might be used 
to separate out or neglect frequency mixing terms 
or harmonic generation. As it stands, no band- 
width restrictions were applied when deriving eqn. 
(|26p - there are only the limitations of the disper- 
sion and/or polarization models to consider. 

5. Mode averaging is where the transverse extent of a 
propagating beam is not explicitly modeled, but is 
subsumed into a description of a transverse mode 
profile; as such it is typically applied to situations 
involving optical fibres or other waveguides. See 
e.g. [g^ for a recent approach, which goes beyond 
a simple addition of a frequency dependence to the 
"effective area" of the mode, and generalizes the 
effective area concept itself. 

A wave equation like that derived above, but limited 
to describing propagation in optical fibres (i.e. a disper- 
sive and third order nonlinear material) , has already been 
studied [s^ ; but it did not consider the effects of diffrac- 
tion or angle dependent refractive index, vector polariza- 
tion terms, or the divergence of P^. It did, however, show 
a stunning level of agreement between uni-directional en- 
velope and PSSD [21| Maxwell equations simulations in 
the case of optical carrier wave shocking - even though 
it described the pulse using an envelope! 
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If desired, we can easily recover wave equations that 
match the SEWA and SVEA wave equations already in 
common use, by applying bandwidth constraints to our 
field, and making approximations based on them. First, 
we set fco = wo(l + S)/c, with S = {uj — ujq)/ujq. Then 
assume that our field has a bandwidth much smaller 
than the carrier frequency ujq, so that E{uJo{l + S)) is only 
non-negligible for (5^1; thus we can now assume S'^ ~ 
0. This bandwidth constraint amounts to an assumption 
about the smoothness of the pulse in the time domain. 
The /cq factor now simplifies to fcg ~ uJq{1 + 2S)/c'^, and 
hence we get a non- envelope but otherwise SEWA-like 
wave equation [T^, which is 



dJ+{Lo) - +z(/3(c.) - kf) E+ioj) + ^^VlE+ioj) 



1 + 2 



UJ — UJQ 



UJQ 



P,{E+{lu),lu). 



(27) 



The next level of bandwidth-limiting approximation 
takes us back to an equation matching the venerable 
SVEA. To achieve this we take such narrow-band fields 
that we can set (5 ~ 0, and so 



dj+iu:) = +1 (/3(w) - kf) E+{lo) + ^^^lE+iLo) 



P,{E+[uj),uj). 



(28) 



Neither of these (SEWA-like, SVEA-like) wave equations 
are required to incorporate an envelope-carrier descrip- 
tion of the fields, or a co-moving frame as demanded by 
the usual SEWA or SVEA derivations; the moving frame 
specified by fc/ above is a mere convenience, and kf may 
be set to zero. Strictly speaking, to match the SEWA 
or SVEA wave equations most closely, we should also set 
/? to a fixed value, and put all of the remaining linear 
dielectric properties of the material into ■ 

Even in the SVEA limit, the factorization technique 
allows us to recover the same propagation equations as 
derived using standard approaches, but this derivation 
now gives us a better (and much simpler) basis on which 
to judge their robustness to strong nonlinearity, angle de- 
pendent refractive indices, and diffraction or transverse 
effects. Note in particular that the linear constraints 
given in section |V] depend only on the material prop- 
erties, and not on the field in any way. The nonlinear 
constraints are the same, but with an additional depen- 
dence on the peak field strength - but importantly, not 
its smoothness or bandwidth. 

It is important to remember that introducing an enve- 
lope and carrier representation of the pulse remains use- 
ful. This is because a well chosen carrier frequency uji 
will almost certainly provide an envelope smoother than 
the field itself; this will provide a more intuitive picture 
but will also have advantages for numerical computation. 



VII. EXAMPLES 

A. Third order nonlinearity 

Third order nonlinearities are common in many mate- 
rials, e.g. in the silica used to make optical fibres ,36i] . 
Here I study propagation in a comparable material, but 
also allow for magnetic dispersion. The propagation is 
based around a wave vector reference /?, where the resid- 
ual frequency dependence of the material refractive index 
is represented by a dimensionless parameter k dependent 
on the linear dispersive parts of the permittivity and 
permeability /^d, so that k = uj{edlJ.dY^^ / P — 1. The in- 
stantaneous electric third order nonlinearity is x^^"*- For 
plane polarized fields, the uni-directional wave equation 
for E^{uj) can be derived from eqn. (P^ . and with the 
usual fco = w/c is 



d,E+ = +ip [1 + k] E+ 



ik^fiL 



2/3 



X^'^Elit)Etit) 



(29) 



where ?"[...] is the Fourier transform that converts the 
time-domain nonlinear polarization into its frequency do- 
main form. 

This is a generalized nonlinear Schrodinger (NLS) 
equation, but is for the full field (i.e. uses no envelope 
description) and retains the full nonlinearity (i.e. retains 
the third harmonic generation term). The only assump- 
tions made are that of transverse fields, weak dispersive 
corrections k, and weakly nonlinear response; these all al- 
low us to decouple the forward and backward wave equa- 
tions. This decoupling then allows us, without any extra 
approximation, to reduce our description to one of for- 
ward only pulse propagation. The specific example cho- 
sen here is for an instantaneous cubic nonlinearity, but it 
is easily generalized to non-instantaneous cases or other 
scalar nonlinearities. 

We can transform eqn. into a NLS equation by 

representing the field in terms of an envelope and carrier, 
where the carrier has a fixed frequency ui and wavevector 
fci; i.e. using 



E+(t) = A{t) exp [i (ujit - kiz)] 

+ A*{t) exp [-t {ujit - kiz)] 



(30) 



In the frequency domain an arbitrary frequency cj dif- 
fers from the carrier frequency wi by an offset A; i.e. 
UJ = uji + A; hence the frequency domain counterpart 
to A{t) is best written ^(A), not A{uj). We proceed by 
setting (3 to have the constant value fci, and ignoring 
the off-resonant THG term, which is usually very poorly 
phase matched. After separating into a pair of complex- 
conjugate equations (one for A, one for A*), this gives 
us the expected NLS equation with diffraction. The cho- 
sen carrier effectively moves us into a frame that freezes 
those carrier oscillations, but this differs from one that 
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is co-moving with the pulse envelope, i.e. one moving at 
the group velocity Vg = doj /dk. After we transform into 
a frame co-moving with the group velocity, where at wi 
we have Kg — u!i{v~^ —v~^), the frequency domain wave 
equation is 



d^A = +tK{A)A 



2k 



x'^^^\A{t)\'A{t) 



2k 



A, 



(31) 



with /^(A) — kn{uji+/S)+Kg. All that has been assumed 
to derive this equation is uni-directional propagation and 
negligible third harmonic generation. This eqn. pip is 
for a magnetically dispersive system broadly comparable 
to that giving rise to the eqn. (12) of Scalora et al. jsHj 
(henceforth eqn. (S12)^); although I have additionally 
retained diffraction and any order of dispersion. 

Many instances of NLS-type equations, such as that of 
eqn. (S12) or simpler forms (e.g. [SGj), are written in the 
time domain, which means that it is more complicated to 
represent the full range of the dispersive response. When 
transforming eqn. (j3ip into the time domain, the dis- 
persion term if (A)A(A) becomes a convolution - but it 
can also be represented as a Taylor series in time deriva- 
tives. This Taylor series is usually reduced to a few low 
order terms, and when using the correct group velocity, 
the lowest order term is a quadratic. Also often seen 
in NLS equations is the self-steepening term (again see 
eqn. (S12)). This self-steepening term be obtained from 
eqn. by expanding k^ = w^/c^ = (wi -I- A)2/c^, 

in a similar manner to deriving a SEWA-like equation 
as discussed in the previous section. Then the leading 
term (oc uj"^) gives the usual nonlinear term, whilst the 
first order contribution (oc 2wiA) gives the single time 
derivative needed for self-steepening in the time domain. 
Also present in eqn. (S12), but not in eqn. (I5T|) is a 
term proportional to x^^-* squared. Here such a term is 
not present because it is second order correction, whilst 
the uni-directional approximation applied here is first or- 
der. Whilst it is possible to incorporate higher-order cor- 
rections, one has to be careful to remain consistent, and 
not miss other significant corrections of the same order, 
nor to include unnecessary terms which should strictly 
be considered negligible. 



B. Second order nonlinearity 

The case of second order nonlinearity is a little more 
complicated, since it typically couples the two possible 
polarization states of the field together [s^ . For simplic- 
ity, I will avoid an exhaustive, detailed derivation from 
first principles, and instead just give example wave equa- 



^ Note that eqn. (S12) has scaled both the time and space param- 
eters. 



tions directly. Indeed, they can be easily inferred directly 
from the format of the coupling in standard treatments. 

In second-order nonlinear interactions such as optical 
parametric amplification (OPA) in lithium borate (LBO) 
using birefringent phase-matching, two field polarizations 
need to be considered. To model the cross-coupling be- 
tween the orthogonally-polarized fields, it is necessary 
to solve for both field polarizations; and to allow for 
the birefringence we need a pair of linear responses, i.e. 

K^{lo), Ky{uj). 

Since it is convenient, I split the vector form of the 
wave equation up into its transverse x and y components. 
The propagation is based around a wave vector reference 
/?, where the residual frequency dependence of the ma- 
terial refractive index in the x ot y directions is repre- 
sented by a dimensionless parameter Ki, for i € {x,y}. 
This K-i is dependent on the linear dispersive parts of 
the permittivities e^.i and permeabilities iid,i, so that 
Ki — u]{ed,i^id,iY^'^ / P — 1- The instantaneous electric sec- 
ond order nonlinear coefficient is x^'^^ ■ Based on eqn. 
(j26p . and for second harmonic generation in the orthogo- 
nal polarization (i.e. a type I OPA), the wave equations 
for E^{uj) and E:^{uj) (with the usual ko = u/c) are 



+ip [1 -h Ei 



^kl^lL 
2/3 



2x'^'^E+{t)Et{t) 



2/3 



-El 



d,E+ = +Z/3 [1 + Ky] E+ 



2/3 



2/3 y 



(32) 
(33) 



where is the Fourier transform that converts the 

time-domain nonlinear polarization into its frequency do- 
main form. The specific example chosen here is easy to 
modify to allow for or incorporate other x*-^-* processes. 
Remarkably, it is also strikingly similar in appearance 
(although not in detail) to the usual SVEA equations 
used to propagate narrow-band pulse envelopes; despite 
the lack of a co-moving frame, and even though they are 
for the field, not an envelope. 

We can transform eqns. ([5^ . ([55]) into the usual equa- 
tions for a parametric amplifier by representing the x and 
y polarized fields in terms of three envelopes and carrier 
pairs: 

Ex{t) = Ai{t) exp [i {ujit — kiz)] 

+ A{ {t) exp [-1 [uoit ~ kiz)] 
+ A2{t) exp [i {uj2t - k2z)] 

+ A*2{t)eyi-p[-i{<jj2t'k2z)] (34) 
Ey{t) = A'iit) exp [i (ust - k^z)] 

+ Al (t) exp [-1 {uj^t - fcgz)] (35) 

where = lui + uj2- After separating into pairs of 
complex-conjugate equations (one each for Ai, one for 
A*), and ignoring the off- resonant polarization terms. 



11 



FCHHG 



Optical pulse propagation 



Dr.Paul.Kinsler@pliysics.org 
http: / / www. kinsler . org /physics / 



Just as for the NLS example above, we also transform 
into a frame co-moving with the group velocity, although 
here we select the group velocity of a preferred frequency 
component (perhaps ws), with e.g. Kg — uj3{v~^ — v~^). 
Choosing (3 for each equation differently, i.e. with /3 e 
{fci, fc2, fcs}, the wave equations for the Ai(uj) are 



2ki 



iK2{A)A2 



2k2 



iK3{A)A3 
ikliiL 



2k. 



2x^^^Asm*2it) 



2/'^A3{t)AUt) 



'/^^A,{t)A2{t) 



^ — lAkz 



(36) 



-lAkz 



2k2 



A2 

(37) 



+iAkz 



2k. 



3- 

(38) 



Here Ki{A) = kiK^{LUi + A) + Kg, with i e {1,2}; and 
K3{A) = k^Ky{ijj^ + A) + Kg. The phase mismatch term 
is Ak — /c3 — fc2 ~ ^1- The only approximations used 
to derive these equations are uni-directional propagation 
and negligible off-resonant polarization terms. 



magnetic responses are allowed. After factorizing the 
second order wave equation into an exact bi-directional 
model, it applies the same slow-evolution approximation 
to all non-trivial effects (e.g. nonlinearity, diffraction), 
and so reduces the propagation equations to a first order 
uni-directional wave equation. My derivation contrasts 
with typical approaches, which often rely on a co-moving 
frame and a sequence of different approximations, such 
as ad- hoc assumption of negligible second derivatives. In 
the appropriate limits, it turns out that many existing 
derivations have given similar but more restricted results 
to those presented here. As a result, with minimal adjust- 
ment, existing numerical and theoretical models could 
be adapted to take advantage of this sounder theoretical 
basis, more straightforward approximations, and simpler 
error-term calculations. 

The improved "factorization" derivation presented 
here allows a term-to-term comparison of the exact bi- 
directional theory with its approximate uni-directional 
counterpart, so that the approximation used (and its 
consequences) is much more easily understood. This 
means that pulse propagation models in the extreme ul- 
trafast and wide-band limits can be made more robust 
- since differences between exact bi-directional and ap- 
proximate uni-directional propagation can be straightfor- 
wardly computed. 
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Appendix A: Factorizing 

Here I present a simple overview of the mathematics 
of the factorization procedure, since full details can be 
found in [sj] . In the calculations below, I transform into 
wave vector space, where the z-derivative dz is converted 
to ik. Also, we have that 0^ = n^iJ^ /<? , and the unspec- 
ified residual term is denoted Q. The second order wave 
equation can then be written 



[dl +P^]E = -Q 
[-e +P^]E = -Q 



E = 



:Q = 



1 



2/3 



1 



(fc-/3) (fc + /3) 
1 



fc + /3 k- P 



Q. 



(Al) 
(A2) 

(A3) 
(A4) 



Now (fc — (3)~^ is a forward-like (Green's function) prop- 
agator for the field, but note that in my terminology, 
it evolves the field. The complementary backward-like 
propagator is {k + (3)~^. As already described in the 
main text, we now write E ~ E^ + E^ , and split the 
two sides up to get 



E+ + E- 
E^ 

[k T /3] E^ 
ikE^ 



2(3 



1 



1 



k + fi k- 13 



Q 



±1 1 



2/3 /fc=F/3 
1 1 



= ± 



2/3fcT/3 



Q 



±^PE^ ± -Q 



(A5) 
(A6) 
(A7) 
(A8) 



Finally, we transform the wave vector space ik terms back 
into normal space to give z derivatives, resulting in the 
final form 



dzE^ 



(A9) 



Appendix B: The no accumulation approximation 

In the main text, I describe the no accumulation ap- 
proximation in spectral terms as a RWA approximation. 
However, it is hard to set a clear, accurate criterion for 
the RWA approximation to be satisfied in the general 
case, since it requires knowledge of the entire propaga- 
tion before it can be justified. In this appendix, I take 



F--^-'^E^ 



(Bl) 



where as noted k can be difficult to determine, and may 
even change dynamically; here we can assume it corre- 
sponds to the propagation wave vector that would be seen 
at if all the conditions holding at a chosen position also 
held everywhere else. On this basis, we can even define 
k = k{z), where by analogy to the linear case we might 
assert that k'^{z) = (3"^ + Q{z)/E{z), so that for small Q, 
we have kE ~ ^{E + Q/2p'^). 

Let us start by assuming our field is propagating and 
evolving forwards (only), with perfectly matched E^ 
fields; so that E^ = £,E^ . but then it happens that Q 
changes by (5Q over a small interval Sz, likewise (, changes 
by S^. The E^ will no longer be matched, and now the 
total field splits into two parts that evolve in opposite di- 
rections. The part that continues to evolve forward has 
E'^ nearly unchanged, but the forward evolving E~ has 
changed size (and is now oc — 6£,)) to stay perfectly 
matched according to the new Q. The rest of the old 
(oc (5^) now propagates backwards, taking with it a tiny 
fraction of the original E^ (and is cx £,S£,). 

Comparing the two backward evolving E^ components 
at z and z + Sz, and taking the limit Sz enables us 
to estimate that the backward evolving E~ field changes 
according to 



dzE, 



2/3 



0. backward 



[k + pY 



IdME, 



0, forward' 



(B2) 



Using the small-Q approximation for fc, we can write 

1 



(fc + /3)- 



[dzQ]i 



(B3) 



where the exponential part removes any oscillations due 
to the linear part of Q; i.e. if Q = xE then 



d E" 

^ G.hackward 



1 



(fc + /3)- 



(B4) 



So here we see that backward evolving fields are only gen- 
erated from forward evolving fields due to changes in the 
underlying conditions (i.e. either material response or 
pulse properties), but that for the reflection to be strong 
those changes will have to be significant on the order of 
a wavelength, or be periodic so that phase matching of 
the the backward wave could occur. 
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